In [91]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
plt.rcParams["animation.html"] = "jshtml"
In [92]:
# Initial conditions and parameters
g = 9.81
l1, l2 = 1.0, 1.0
m1, m2 = 1.0, 1.0
dt = 0.01
T_max = 1000
N = int(T_max / dt)
t = np.linspace(0, T_max, N)
y0 = np.radians([120, -10, 0, 0])
In [93]:
# Functions for dynamics
def alpha_1(theta_1, theta_2):
return (l2/l1) * (m2/(m1 + m2)) * np.cos(theta_1 - theta_2)
def alpha_2(theta_1, theta_2):
return (l1/l2) * np.cos(theta_1 - theta_2)
def f1(theta_1, theta_2, w2):
return -(l2/l1) * (m2/(m1 + m2)) * (w2 ** 2) * np.sin(theta_1 - theta_2) - (g/l1) * np.sin(theta_1)
def f2(theta_1, theta_2, w1):
return (l1/l2) * (w1 ** 2) * np.sin(theta_1 - theta_2) - (g/l2) * np.sin(theta_2)
def RK_term(t, y, dt, f):
k1 = dt * f(y, t)
k2 = dt * f(y + 0.5 * k1, t + 0.5 * dt)
k3 = dt * f(y + 0.5 * k2, t + 0.5 * dt)
k4 = dt * f(y + k3, t + dt)
return (1/6) * (k1 + 2 * k2 + 2 * k3 + k4)
def domega(y, t):
theta_1, theta_2, w1, w2 = y
c = 1 / (1 - alpha_1(theta_1, theta_2) * alpha_2(theta_1, theta_2))
g1 = c * (f1(theta_1, theta_2, w2) - alpha_1(theta_1, theta_2) * f2(theta_1, theta_2, w1))
g2 = c * (f2(theta_1, theta_2, w1) - alpha_2(theta_1, theta_2) * f1(theta_1, theta_2, w2))
return np.array([w1, w2, g1, g2])
def T_energy(point):
theta_1, theta_2, w1, w2 = point
T = 0.5 * (m1 + m2) * (l1 ** 2) * (w1 ** 2) + 0.5 * m2 * (l2 ** 2) * (w2 ** 2) + m2 * l1 * l2 * w1 * w2 * np.cos(theta_1 - theta_2)
V = -(m1 + m2) * g * l1 * np.cos(theta_1) - m2 * g * l2 * np.cos(theta_2)
return T + V
def RK_4_DP():
y = [y0]
for i in range(1, N):
y.append(y[i - 1] + RK_term(t[i - 1], y[i - 1], dt, domega))
total_energy_rk = np.array([T_energy(p) for p in y])
return y, total_energy_rk
def Euler_DP():
y = [y0]
for i in range(1, N):
y.append(y[i - 1] + dt * domega(y[i - 1], t[i - 1]))
total_energy_eu = np.array([T_energy(p) for p in y])
return y, total_energy_eu
In [94]:
rk_y, rk_E = RK_4_DP()
eu_y, eu_E = Euler_DP()
In [95]:
# Energy comparison
plt.figure(figsize=(10, 5))
plt.plot(t, rk_E, label="RK4 Energy", color='green')
plt.plot(t, eu_E, label="Euler Energy", color='orange')
plt.axhline(rk_E[0], linestyle='--', color='red', alpha=0.5, label="Initial Energy")
plt.title("Energy Conservation: Runge-Kutta vs Euler")
plt.xlabel("Time (s)")
plt.ylabel("Total Energy (J)")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
In [96]:
# Phase portraits
rk_theta_1 = np.array([s[0] for s in rk_y])
rk_theta_2 = np.array([s[1] for s in rk_y])
eu_theta_1 = np.array([s[0] for s in eu_y])
eu_theta_2 = np.array([s[1] for s in eu_y])
plt.figure(figsize=(10, 5))
plt.plot(rk_theta_1, rk_theta_2, label="RK4", color='green')
plt.plot(eu_theta_1, eu_theta_2, label="Euler", color='orange', alpha=0.6)
plt.title(r'$theta_1$ vs $theta_2$ (Phase Portrait)')
plt.xlabel(r'$theta_1$ (rad)')
plt.ylabel(r'$theta_2$ (rad)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
In [97]:
rk_omega_1 = np.array([state[2] for state in rk_y])
rk_omega_2 = np.array([state[3] for state in rk_y])
eu_omega_1 = np.array([state[2] for state in eu_y])
eu_omega_2 = np.array([state[3] for state in eu_y])
plt.figure(figsize=(8,6))
plt.plot(rk_omega_1, rk_omega_2, label="RK4 ω", color='green')
plt.plot(eu_omega_1, eu_omega_2, label="Euler ω", color='orange', alpha=0.6)
plt.title(r'$\omega_1$ vs $\omega_2$ (Angular Phase Portrait)')
plt.xlabel(r'$\omega_1$ (rad/s)')
plt.ylabel(r'$\omega_2$ (rad/s)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
In [98]:
# Animation functions
def get_double_pendulum_coords(theta1, theta2):
x1 = l1 * np.sin(theta1)
y1 = -l1 * np.cos(theta1)
x2 = x1 + l2 * np.sin(theta2)
y2 = y1 - l2 * np.cos(theta2)
return x1, y1, x2, y2
fig, ax = plt.subplots()
x1_0, y1_0, x2_0, y2_0 = get_double_pendulum_coords(rk_theta_1[0], rk_theta_2[0])
line, = ax.plot([0, x1_0, x2_0], [0, y1_0, y2_0], lw=2, c='k')
bob1 = ax.add_patch(plt.Circle((x1_0, y1_0), 0.05, fc='r'))
bob2 = ax.add_patch(plt.Circle((x2_0, y2_0), 0.05, fc='b'))
ax.set_xlim(-1.5*(l1+l2), 1.5*(l1+l2))
ax.set_ylim(-1.5*(l1+l2), 1.5*(l1+l2))
ax.set_aspect('equal')
ax.grid(True)
In [99]:
def animate(i):
x1, y1, x2, y2 = get_double_pendulum_coords(rk_theta_1[i], rk_theta_2[i])
line.set_data([0, x1, x2], [0, y1, y2])
bob1.set_center((x1, y1))
bob2.set_center((x2, y2))
return line, bob1, bob2
ani = animation.FuncAnimation(fig, animate, frames=len(rk_theta_1), interval=dt * 1000)
%matplotlib notebook
plt.rcParams["animation.html"] = "jshtml"
ani
Animation size has reached 20981281 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.
Out[99]: